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Abstract 

We present a novel computational methodology for solving the scalar nonlinear Helmholtz equation 
(NLH) that governs the propagation of laser light in Kerr dielectrics. The methodology addresses 
two well-known challenges in nonlinear optics: Singular behavior of solutions when the scattering in 
the medium is assumed predominantly forward (paraxial regime) , and the presence of discontinuities 
in the optical properties of the medium. Specifically, we consider a slab of nonlinear material which 
may be grated in the direction of propagation and which is immersed in a linear medium as a 
whole. The key components of the methodology are a semi-compact high-order finite-difference 
scheme that maintains accuracy across the discontinuities and enables sub-wavelength resolution 
on large domains at a tolerable cost, a nonlocal two-way artificial boundary condition (ABC) 
that simultaneously facilitates the reflectionless propagation of the outgoing waves and forward 
propagation of the given incoming waves, and a nonlinear solver based on Newton's method. 

The proposed methodology combines and substantially extends the capabilities of our previous 
techniques built for ID and for multi-D. It facilitates a direct numerical study of nonparaxial 
propagation and goes well beyond the approaches in the literature based on the "augmented" 
paraxial models. In particular, it provides the first ever evidence that the singularity of the solution 
indeed disappears in the scalar NLH model that includes the nonparaxial effects. It also enables 
simulation of the wavelength-width spatial solitons, as well as of the counter-propagating solitons. 
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1. Introduction 



1.1. Mathematical Models 



The propagation of electromagnetic waves in materials is governed by Maxwell's equations with 
appropriately chosen material responses. The responses characterize the dependence of material 
properties — magnetic permeability, electric permittivity, and conductivity — on the location and 
frequency of the propagating field. For high intensity radiation, the material quantities may also 
depend on the magnitude of the propagating field, which makes the responses nonlinear. 

In nonlinear optics, one is often interested in studying the propagation of monochromatic waves 
(continuous-wave laser beams) through transparent dielectrics. In this case, the generation of 
higher harmonics and nonlinear coupling between different (temporal) frequencies can often be 
neglected, and accordingly, a time-harmonic solution can be assumed. The magnetic field can 
then be eliminated, and Maxwell's equations transform to a second-order differential equation with 
respect to the electric field, known as the vector Helmholtz equation, see [H. If the material is 
isotropic and, in addition, the electric field is assumed linearly polarized, then one arrives at the 
scalar nonlinear Helmholtz equation (NLH): 



to, 



0, n'^{x, \E\) = nl{x) + 2nQ{x)n2{x)\E 



2a 



(1) 



where o" > and n is the refraction index. In physical materials one always has a = 1, so that the 
dependence of on \E\ is quadratic. In equation ([T|), x = [xi, . . . , xd] are the spatial coordinates, 
E = E{x) denotes the scalar electric field, loq is the laser frequency, c is the speed of light in 
vacuum, A = d'^^ + . . . + d'^^ is the L'-dimensional Laplacian, no is the linear index of refraction, 
and 712 is the Kerr coefficient. Both no and n2 are assumed real, so that the medium is transparent 
or lossless. The coordinate xd will also be denoted by z and will hereafter be referred to as 
longitudinal, whereas the remaining direction(s) x± = [xi, . . . ,X£)-i] will be called transverse. 

Our primary physical setup involves a slab of Kerr material surrounded on both sides by the 
linear homogeneous medium in which no 



n, 



ext 



and n2 = 0, see Figure 1(a) 
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(a) The three-layer physical setup. 
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(b) The multi-layer physical setup. 



We introduce the linear wavenumber ko 

^ext g^^j _ /(^cxt\2 



u;o7^o^*/c and the normalized quantities I'i^x) 



noix)/nl'' 



2nQ{x)n2{x) / {n'^^) , and then recast equation ([T]) as 
AE{x) + kl 



[v\x) + e{x)\Er] 



E = {). 



(2) 



Note that the Kerr coefficient n2(a;) is always discontinuous at the interface planes z = and 



-^max) see Figure 1(a) The linear index of refraction no(a;) may also be discontinuous at the 



interface planes. Discontinuities in no(a;) and n2{x) immediately give rise to those in i^{x) and e{x) 
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see equation ( 
Kerr material 



Thus, for the typical experimental setting that involves a slab of homogeneous 



the coefficients of equation ([2]) are piecewise-constant: 



1, 

h 



z <0, 
0< z< Z„ 



and 



e{z, x±] 



0, z <0, 

e™*, < Z < ^max) 
0, Z > -^max- 



(3) 



Discontinuities in the coefficients ([3]) imply that additional conditions will be required for the 
NLH ([2]) at the interfaces z = and z = ^max- These conditions can be obtained by analyzing the 
corresponding Maxwell's equations. They reduce to the continuity of the field E{z) and its first 
normal derivative see Appendix |Al When building a numerical approximation, the presence of 
material discontinuities requires special attention (Sections [5] and ED . 

The problem is driven by a laser beam that impinges on the Kerr material from the outside and 



causes a local increase in the overall index of refraction as it propagates through, see Figure l(a 



Since light rays bend toward the areas with higher refraction index, the impinging beam self-focuses 
inside the Kerr medium. The material discontinuities at 2: = and z = Z^^ax reflect a portion of 
the forward propagating wave, resulting in backward propagating waves. Moreover, the nonlinearly 
induced nonuniformities of the refraction index may also scatter the radiation backwards. The 
presence of waves propagating in opposite directions implies that the boundary conditions for the 
NLH ([2]) must ensure the reflectionless propagation of all the outgoing waves (regardless of their 
direction of travel and the angle of incidence at the outer boundary) and at the same time correctly 
prescribe the given incoming beam at the boundary, see Figure [H Such boundary conditions are 
called two-way boundary conditions i, see Sect.onifor ID and Section^for multi-D. 
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Figure 1: Schematic of the boundary conditions in the longitudinal direction: A) One-way radiation boundary 
condition at z = Zy^ax + S; B) Two-way radiation boundary condition at z = —S. 



One can also consider a simplified model that would account only for the forward propagating 
component of the field. Let z = he the direction of the impinging laser beam, and let us 
also consider the simplest case of = 1 and e = e™* inside the Kerr medium. Then, introducing 
the ansatz E = e^^°^(p, where (p = <^(^) is assumed to vary slowly compared with the fast carrier 
oscillation e^^°^, one can neglect the small (pzz term (paraxial approximation), and reduce the 
NLH dJ]) to the nonlinear Schrodinger equation (NLS): 

2iko(l)z{z,x^)+A^(l) + kie\(l)\^''(l) = 0, 0<z<Z,^ax, (4) 



^ This setup withstands an easy generalization to the case of multiple plane-parallel layers, see Section [1.31 
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which governs the envelope (f). In equation A_l = + . . ■ + d'^^_^ denotes the transverse Lapla- 
cian. The NLS ^ supports only forward propagation because the assumption of slow variation of 
does not leave room for any ~ e~^^^^ components in the solution. Equation @ is first order in 
z and, unlike the NLH, requires a Cauchy problem to be formulated and solved with the "initial" 
data provided by the impinging wave and specified, say, at z = (see, e.g., 0, 0] for detail). 

It is well known that solutions of the NLS (j4]) exist globally when a{D — 1) < 2, the sub- 
critical NLS, but can become singular, i.e., collapse at finite propagation distances, when either 
a{D — 1) > 2, the supercritical NLS, or a (L» - 1) = 2, the critical NLS [1]. As shown by Wein- 
stein 0], a necessary condition for singularity formation in the critical NLS is that the input power 
exceeds the critical power P^. The value of Pc is equal to the power of the ground-state solitary 
wave solution of the NLS; it can be calculated analytically for D = 2 and numerically for D > 2. 

A question that has been open in the literature for over forty years is whether the more compre- 
hensive NLH model for nonlinear self-focusing eliminates the singular behavior that characterizes 
collapsing solutions of the critical and supercritical NLS. Unfortunately, the fundamental issue of 
solvability of the NLH and regularity of its solutions still remains unaddressed for many important 
settings. Only the one-dimensional case, when equation ([2]) becomes an ODE, has been studied 
extensively^, and exact solutions have been obtained using a combination of analytical and numerical 
1 S U; 12, 13]. In multi-D, there have been indications that solutions of the NLH 



means 



may exist even when the correspon ding NLS solutions become singular, based on both numerical 



study of "modified" NLS equations [IJ, lla, |l6(] , and on asymptotic analysis [17|, but these studies 
did not account for backscattering. Recently, Sever employed a Palais-Smale type argument and 
has shown that the multi-D NLH is solvable in the sense of and that the solution is not 
unique [3]. His argument, however, only applies to self-adjoint operators, whereas the physical 
setups considered in this study require radiation boundary conditions. 

1.2. Numerical Method 

The new coinputational methodology for the NLH that we present builds up on our previous 
work 0, Q, 21] and extends it substantially. We introduce a new semi-compact discretization 



and a new Newton's solver, and the ensuing capabilities include an explicit demonstration of the 
removal of singularity that "plagues" the NLS, and the computation of narrow nonparaxial solitons. 

Specifically, we solve the NLH ([2]) for two different cases. The first one corresponds to the 
critical NLS {cr{D — 1) = 2). We consider both the two-dimensional quintic nonlinearity D = 2 and 
(7 = 2 (planar waveguides), and the three-dimensional cubic nonlinearity D = 3 and a = 1 (bulk 
Kerr medium, for which we additionally assume cylindrical symmetry). As a{D — 1) = 2 for either 
setting, one can expect that the role of nonparaxiality and backscattering will be similar. This study 
goes beyond the investigation of the "modified" NLS's |li,[li,[li, Il7| ]. and the results reported 



m 



Section 17.21 provide the first ever numerical evidence that the collapse of focusing nonlinear waves 
is indeed arrested in the NLH model, which incorporates the nonparaxiality and backscattering. 

The second case we analyze is that of a planar waveguide with cubic nonlinearity (D = 2 and 
(7 = 1). In this subcritical case, solutions to the NLS do not collapse. Instead, the laser beam can 
propagate in the Kerr medium over very long distances without changing its profil^ — the type 
of behavior often referred to as spatial soliton. Solitons have been studied extensively as solutions 
to the NLS. For beams that are much wider than the optical wavelength, it is generally expected 



^ In this case, self-focusing balances diffraction exactly. 



4 



that the "subcritical" NLH will have similarly looking solutions. However, it was not until our 
paper 20(] that it has become actually possible to study the effect of nonparaxiality and backscat- 
tering on solitons. The methodology proposed in this paper allows us to go further and demonstrate 
numerically the existence and sustainability over long distances of very narrow spatial solitons for 
the NLH, basically as narrow as one carrier wavelength A = 27r/fcoi see Section [7.11 Furthermore, 
the NLH appears particularly well suited for modeling interactions between counter-propagating 
solitons, as a boundary value problem can naturally be formulated. In the NLS framework, on the 
other hand, the two counter-propagating solitons will imply two opposite directions of marchingH 
The discrete approximation of the NLH must be high order so as to minimize the number of 
points per wavelength required for solving equation ([2]) with sub-wavelength resolution on a large 
domain, and for resolving the small-scale phenomenon of backscattering against a background of the 
forward-propagating wave. It must also maintain its accuracy across the material discontinuities. 
As the geometry is simple, and the discontinuities are only in the longitudinal direction, we can 
approximate the NLH by finite differences on a rectangular grid. In the case D = 2, it will be a 
Cartesian grid of coordinates (x, z). In the case D = 3, we still want to have only two independent 
spatial variables and hence employ cylindrical symmetry. The NLH ([2]) is then approximated on 
the rectangular grid of cylindrical coordinates {p,z), where p = (x^ + y^)^/^. In doing so, the 
discontinuities that are confined to transverse planes will always be aligned with the grid. 



B 




Figure 2: Stencils in 2D: A) Standard central difference fourth-order stencil, as in our previous work 0, [l^, HJ]; B) 
Compact 3x3 fourth-order stencil for linear operators, as, e.g., in [2^, C) Semi-compact stencil used in this 
work. 



In our work [2, l20[] , we used the standard fourth-order central differences (five node stencil in 
each coordinate direction) to approximate the NLH ([2|) on a rectangular grid, see Figure [2]A. While 
this approach works well in the regions of smoothness, it deteriorates to second-order accuracy in 
regions of material discontinuities. In the recent paper [l^, we discretized the one-dimensional 
NLH with fourth-order accuracy using compact finite volumes and a three node stencil. This 
discretization handled the material discontinuities with no deterioration of accuracy and was also 



^ Counter-propagating beams have been simulated using two coupled NLSs [221, but this approach involves some 
approximations which are not needed in the NLH, and whose validity is unclear. 
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extended to higher orders in the Unear case [25|]. However, the extension of the scheme of [191 ] to 
multi-D is not straightforward. Therefore, in the current paper we adopt a hybrid approach. We use 
the standard fourth-order central differences in the transverse direction, and a compact fourth-order 
finite difference discretization on three nodes in the longitudinal direction, see Figure ElC. 

The five node transverse stencil does not impair the accuracy because there are no discontinuities 
in that direction. The three node longitudinal part of the scheme is supplemented by one-sided 
differences that implement the required interface conditions at the points of discontinuity. In doing 
so, the compact stencil eliminates the need to use those special differences anywhere except at 
the discontinuities themselves. Another advantage of having a three node compact stencil in the 
longitudinal direction is that it leads to matrices with a narrower bandwidth. 

The interior discretization is supplemented by nonlocal two-way artificial boundary conditions 
(ABCs) set at z = —5 and z = Z^ax + see Figure [H and by local radiation boundary conditions 



at the transverse far-field boundaries. The discrete ABCs are similar to those of [20|], but having 
a three node compact stencil greatly simplifies their construction because, unlike in the case of a 
five node stencil, there are no additional evanescent modes in the discretization, see Section 13. 4[ 

The solver employed in [2, [i^, 21] was of a fixed-point type. On the outer iteration loop, the 
nonlinearity in equation ([2]) was frozen, and a linear Helmholtz equation with variable coefficients 
was obtained. This linear equation was then solved iteratively on the inner loop, essentially by 



building a sequence of Born approximations [26|] . This double- loop iterative method was shown to 
converge for (subcritical) solitons and for "mild" critical cases, but has never been able to produce 
convergent solutions for incoming beams that become singular in the NLS model. 

In [3], we have demonstrated that the iterations' convergence in [3, 21] breaks down far 
below the power threshold for non-uniqueness of the one-dimensional problem. This suggested that 
the convergence difficulties in 0, 23, 21] were not related to the loss of uniqueness by the solution 



3,8, but rather to the deficiencies of the iteration scheme itself. The latter may be (partially) 
accounted for by the known conver genc e limitations of the Born approximations, because they can 



be interpreted as a Neumann series [27[ for the corresponding integral operator [26l . Section 13.1.4]. 

An alternative iteration proposed in [l^ is based on Newton's method. As, however, the Kerr 
nonlinearity is Frechet nondifferentiable for complex-valued E, for Newton's method to apply the 
NLH has to be recast as a system of two equations with real unknowns. The one-dimensional 



numerical experiments of [l9j] demonstrate robust convergence of Newton's iterations for a wide 
range of input powers. Therefore, in this paper we implement Newton's method for solving the 
multi-dimensional NLH ([2]), see Sectional As shown in Section [721 the method converges for initial 
conditions that lead to singularity formation in the critical NLS model, for both D = 2 and D = 3. 

1.3. Extension to the Multi-Layer Case 



Instead of having a homogeneous Kerr material in the nonlinear region as shown in Figure l(a 



we can analyze the case of a layered (grated) material as shown in Figure 1(b) In doing so, the 
linear material outside of the Kerr slab still remain homogeneous. 

The corresponding extension of the mathematical model is straightforward. It amounts to 
introducing a fixed partition of the interval [0, .^max]: 

Q = Zq < ■ ■ ■ < Zi < ■ ■ ■ < ZL = ^max, (5a) 

so that the material characteristics are constant within each sub-interval: 

v{z,Xi_) = ui, e{z,x^) = ei, for ze{zi,zi+i), (5b) 
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whereas at the interfaces (jSap they may undergo jumps. Altogether, this leaves the coefficients of 
equation ([2]) piecewise constant in z. The additional interface conditions required by equation ([2]) 
are the same as before — continuity of E and ^ between the layers, see Appendix Rl 

1. J^. Structure of the Paper 

In Section [51 we illustrate the main concepts of the continuous formulation and the discretization 
for the one-dimensional NLH. In Section [3l we describe the continuous formulation of the problem 
and the discretization for the two-dimensional Cartesian NLH and for the three dimensional NLH 
with cylindrical symmetry. In Section U we introduce Newton's solver for the resulting system of 
nonlinear equations on the grid. Section [5] provides a summary on the numerical method. Section [6] 
relates the input beams for the NLH and the corresponding NLS models, and Section [7] contains 
the results of simulations. Finally, Section [8] presents our conclusions and outlines directions for 
future work. Note also that some of the results shown hereafter were previously reported in 28l |. 
That paper, however, did not contain any description of the numerical method. 

2. The NLH in One Space Dimension 

In this section we consider the one-dimensional NLH with constant material coefficients z/^ and 
e for < 2: < 

Zmaxi which means that there are two discontinuities at z — and z — -^max; 



but no discontinuities in the interior of the Kerr slab, see Figure 1(a) In Section [2.11 we present 
the continuous formulation of the problem, and in Section 12.21 we introduce a compact discrete 
approximation. In Section 12. 3[ we briefly discuss the extension to the multi-layer case outlined is 
Section 11.31 This simple one-dimensional case illustrates the key ideas and notations that will be 
used later in the more complex multi-dimensional cases. 

2.1. Continuous Formulation 

Consider a homogeneous slab of the Kerr material immersed in an infinite linear medium. The 
propagation of the electric field is governed by the ID NLH equation inside the Kerr material: 

"^'^^'^ + ki + 6 \E\''') E = 0, 0<z< Z^ax, (6a) 



dz^ 

and by the linear Helmholtz equation outside the Kerr material: 

+kiE = 0, Z<0 01 Z> Z,-nax. (6b) 

At the material interfaces z = and z = Z^^ax, the field and its first derivative must be continu- 
ous Hi]: 

. s . ^ dE , , dE , , 
E(„+) = EiO-), -(0+) = -(0-), ^^^^ 

dE , s dE 

E{Z-[na.x~\~) = E[Znjax~), '~, — (■^max + ) = ( -^max ~ ) • 

dz dz 

We consider the case of two incoming waves with known characteristics that travel toward the 
Kerr region [0, .^max] from z = — 00 to the right and from z = +00 to the leftl^l The overall field 



® Hereafter, we slightly generalize the schematic depicted in Figure [T] in that we allow for incoming waves to 
impinge on both interfaces, at 2; = and z — ^max- 
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may also have scattered components, which are outg oing with respect to the doma-in [0, ^max 

1 and 

which are not known ahead of time. The general solution to equation fl6bl) outside [0, -^max] ^ su- 
perposition of the right-propagating wave e*'^"^ and the left-propagating wave e"*'^''^. Consequently, 
the field outside [0, Z^^y\ shall be sought for in the form: 



E{z) 



— CX) < -2 < 0, 
^max < Z < OO, 



(7) 



where E^^^ is a given amplitude of the incoming wave that travels to the right from z = —oo 
and impinges on the Kerr medium at z = 0, whereas C\ is the amplitude of the outgoing wave 
traveling to the left toward z = — oo, which is not known ahead of timelll. Likewise, E^^'^'' is a given 
amplitude of the incoming wave that travels to the left from z = +oo and impinges on the Kerr 
medium at z = .^maxj whereas C2 is the amplitude of the outgoing right-traveling wave, which is 
not known ahead of time. 

Representation ([7]) is to be enforced by the ABCs that should prescribe the given values of Ef^^ 
and -E'in™'"' ^^'^ the same time allow for the arbitrary values of Ci and C2. In 19], we have 
set such ABCs precisely at the material interfaces, and have shown that they were given by the 
inhomogeneous Sommerfeld type relations: 



-^ + iko\E 
az 



2ikoE9 



z=0 



^-iko]E 

az 



In this paper, we set equivalent ABCs at a certain distance 5 > away from the interfaces, see 
Figure [U and inside the linear regions: 



^ + iko]E 
az 



dz 



iko 1 E 



2ikoe-"''^^Ei^, 



2ikne~'''°^ eFj^-^ . 



(8) 



2 — '^max+^ 



As we shall see, the separation between the material interfaces z = and z = Zmax and artificial 
boundaries z = —6 and z — -^max ~l~ ^ simplifies the discretization of the problem, because the 
continuity conditions ([6c|) and the boundary conditions ([8|) can be discretized independently of 
each other, see Sections 12.2.21 and 12.2.31 respectively. 



2.2. Discrete Approximation 

The one-dimensional problem ([6]), ([8]) will be approximated using compact fourth-order finite 
differences. We first discuss the discrete approximation of equations (j6ap and (I6b|) . then the ap- 
proximation of the interface condition ([5c|) . and finally the approximation of the two-way ABCs ([5]). 
In what follows, we introduce some notations that will be particularly helpful in multi-D. 



^ Physically, the left-traveling outgoing wave Cie"''""^ may have two sources: A portion of the right-traveling wave 
iSfnt-e**^"^ may get scattered to the left by the Kerr material slab, and a portion of the left-traveling wave i^.^^a^e"'*^"^ 
may be transmitted through by the Kerr material slab. In the nonlinear problem, these phenomena are coupled and 
cannot be easily distinguished from one another. 
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We begin with setting up a uniform grid of + 7 nodes on [—5, ^max + S]'- 

zn = n-h, h=^^, n = -3,-2,...,N + 2,N + 3, (9) 

so that 

Zq = 0, ZN = ^max, ^ = 3/l. 

We also denote by En and P„ = \En\'^"En the values of E and of \E\'^"E at the grid nodes z^- 
Finally, we introduce a special notation D for central difference operators, with the order of accuracy 
in the superscript and the differentiation variables in the subscript. For example, 

def En+l — 2En + En~l d?E 



d^?)e 



+ o[h^). 



/l2 

2.2.1. Approximation of the Equation 

Inside the Kerr medium, i.e., for n = 1,...,A^ — 1, the material coefficients v"^ and e are 
constant, and hence the field E{z) is smooth. Using Taylor's expansion of the field, we obtain from 
the standard second-order central difference approximation: 

Df)E = - 2En + g.-i ^ ^^^^^ ^ ^d.,,,En + O {h') . (10) 

Then, recasting the one-dimensional NLH (j6ap as 

d,,En = (zy2 + e|^„|2-) En = -kl{v^ En + eP„), 

we can approximate the term dzzzzEn on the right-hand side of (llOp with second-order accuracy as 

dzzzzEn = Df^dzzEn + O {h^) = -klof) {u^En + cPn) + O [h^] . 

This yields a compact fourth-order approximation for the second derivative: 

dzzEn = Df)En + ^Df) {u'En + eP„) + O {h'). 
Then, the resulting scheme for the one-dimensional NLH (|6ap at the interior nodes reads: 

Df)En + kl {l + ^^2)) {u'^En + ePn) = 0, ^^^^ 

n= l,...,iV- 1. 



This approach is sometimes called an equation-based approximation [24]. 

Outside the Kerr medium, i.e., for n < and n > N, the foregoing derivation is repeated with 
1/2 = 1 and e = 0, which yields a compact fourth-order approximation of the linear Helmholtz 
equation ([6b]): 

+ ^) Di^jEn + k^,En = 0, ^^^^ 

n = -3,-2,-1, and n = iV + 1, + 2, + 3. 

Note that equation (jl2p for the outermost grid nodes n = —3 and n = iV + 3 will involve the 
ghost values S_4 and Epf^4^, respectively. These ghost values will be determined from the discrete 
two-way ABCs, see Section [2.2.31 
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2.2.2. Approximation at the Interfaces 

At the material interfaces z = and z = ^max (i-e., grid nodes n = and n = N) the discretized 
field is given by £^0 and Ej^, respectively. Hence, the continuity of E at the interface is automatically 
guaranteed, and only the continuity of Ez, see formula (f6c]l . requires special attention. The latter 
is enforced by approximating the derivatives at the interfaces with fourth-order one-sided finite 
differences. We again use the differential equation to eliminate one grid point from the one-sided 
stencil and reduce it from the conventional five nodes to four. While reducing the size of the stencil 
at the interface is not as important as in the interior and exterior of the Kerr material, numerical 
observations show that in some cases it may bring down the truncation error at the interface by a 
factor of two. 

Using Taylor's expansion and the one-dimensional NLH (|6ap . we can write: 



dE 

dz 



-85^0 + 108£;i - 27E2 + ^E^ 3h d'^E 



^^0+ 66/1 11 dz2 

U2! 



z=0+ 



_ -85Eo + W8E^-27E2 + 4E3 3^h , , ^ p \.n(h^\ 

Repeating the calculation for £'2(0—) and equating the resulting approximations for £'2(0—) and 
Ez{^+), we have: 

4S_3 - 27E^2 + 108£;_i - 170£;o + IO8S1 - 27£;2 + ^E^ 

mh 

Qhkl ( vl^ + z^o+ ^ , + eo^ 



+^ [f^^E^ + -^—f^P^) = 0- (13) 

Then, substituting i/q- = 1) £0- = 0, !/o+ = ^ a-iid eo+ = we obtain: 

4£_3 - 27£;_2 + 108£;_i - 170£o + 108£i - 27£;2 + 4£^3 

mh 

Qhkl fl + u"^ e \ , , 

+-rr -^^Eo + -Po)= 0. (14a) 
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A similar equation is obtained for the interface at n = A^: 

4:Ei^_^ — 27Ej\f^2 + 108£jv-i ~ 170£^ + 108£'7v+i ~ 27£jv+2 + ^£^^3 

mh 

+^ [-^En + ^EnJ = 0. (14b) 

2.2.3. Two- Way Boundary Conditions 

At the exterior nodes n < and n > N, the one-dimensional Helmholtz equation ()6bp is approxi- 
mated with fourth-order accuracy by the constant coefficient homogeneous difference equation (I12p . 
This equation can be recast as 

+ k =0, where k = ^^^^^K- (15) 

The general solution of equation (fT5]) is of the form £'„ = C+g" -|- C-q~'^, where 



-|- iyl — r^ and q = r — iyV 



r^ 



10 



are roots of the corresponding characteristic equation — 2/r + q = and r = (1 — k'^h? /2)~^ . 
These roots are complex conjugate and have unit magnitudes. Moreover, they satisfy q = 
^ikoh ^-^ Q (/jS^^ g^^^ q-1 _ Q-tkoh (^\j^Q (/i^)). Hence, the discrete solution approximates 

the right-going wave e^^^'^^ = e**^°^, and the discrete solution q~"' approximates the left-going wave 
^-ikonh _ g-ifco^^ with fourth-order accuracy. Consequently, the discrete counterpart of equation ([7|) 
is 

E = / ^inc9" + Ciq--, -oo < n < 0, 

" \c2g"-^ + Sifr''e"^""^\ iV<n<oo. 

Applying equation (jl6p at n = —3 and n = —4, we can eliminate the unknown constant Ci and 
express the value of the field at the ghost node n = —4 as 

£;„4 = (q-^ - q)q~^Ef^c + 9^-3- (17a) 
Likewise, applying equation (jl6p at n = -|- 3 and n = -|- 4, we obtain: 

En+a = {q~^ - q)q~''E^-- + qE^+s- (17b) 

Relations (jl7p provide a fourth-order accurate approximation to the boundary conditions ([8]) for 
5 = 3h. Relation (jlTah is substituted into equation (fT2|) forn = —3 and relation (|17bp is substituted 
into equation (jl2p for n = A^-|-3. This eliminates the ghost values from scheme (jl2p and closes the 
system of difference equations on the grid ([9]). 

2.3. Extension to the Multi-Layer Material 

In the case of a grated Kerr material described in Section [1.31 there are additional discontinuity 
points defined by formula ()5ap . The interface conditions at each discontinuity point z are the same 
as at z = and z = Z^^x- 

„, X ^ / s dE , . dE , , 
= -(.+) = -(.-). 

Hence, in the simple case when z happens to be at one of the grid nodes, the discrete continuity 
condition at z is given by the same expression as (|13p . If the discontinuity point does not coincide 
with any grid node, one can construct a separate uniform grid for each sub-interval, and the 
extension to the multi-layer case will then be straightforward. 

3. The NLH in Two and Three Space Dimensions 

3.1. Continuous Formulation 

Here, we build a continuous formulation for the case of a homogeneous slab of the Kerr material 



which occupies the region < z < .^maxj see Figure 1(a) As in the one-dimensional setting, we 
will later generalize the method to the multi-layer case, see Section [3.51 

We first consider the two-dimensional Cartesian geometry case x = (z, x). This case models the 
physical case of propagation in planar waveguides, where the dynamics in y can be neglected. In 
this case, the computational domain is truncated in the transverse direction to x £ [— Amaxi Amax]- 
In the longitudinal direction, we truncate the computational domain at a certain distance 6 from 
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the interfaces, to z £ [—5, ^max + 5]. As before, the electric field is governed by the scalar NLH 
equation inside the Kerr medium [cf. equation ([2])]: 

E,,{z,x)+E^^ + kUu^ + e\E\^'')E = 0, 

^ ' (18a) 

(z^x) G (0, ^niax) ^ [ ^max 1 ^max] ) 

and by the linear Helmholtz equation outside the Kerr medium (where v =\ and e = 0): 

E,,{z,x)^E^^^klE = ^, ^^^^^ 

(z, X) G {[-(5, 0) U (^max, ^max + (5]} X [-Xmax, -'^max]- 

At the material interfaces z = and z = Z^^i^^, the field E and its normal derivative Ez are 

continuous for all X £ [— -'^maxi -'^max]: 

E{0+,x)=E{0-,x), E,{0+,x) = E,{0-,x), ^^^^^ 

^(^max+i^;) = E[Zj:ns^^ — ,x), i?^ (^max+j 3;) = i?^ ( Zjnax ~ ) 2; ) . 

We also consider the case of three spatial dimensions, which models the propagation in bulk 
medium. In order to reduce the computational costs, we assume that the field is cylindrically 
symmetric E{x) = E{z,p), where p = \x±\ = \J x^ ^ y^. This enables us to solve the problem with 
only two independent spatial variables. In this case, the computational domain in the transverse 
direction is p G [0, Pmax]) and the scalar NLH equation inside the Kerr medium is 

£;,,(z,p)+^pp + i^p + /c2(i/2 + e|^|2-W = o, 

P ^ ' (19a) 

(z, p) £ (0, Zm^^) X [0 
The linear Helmholtz outside the Kerr medium is 

E,,{z, p) + Epp + -Ep + klE = 0, 

P (19b) 

{z,p) G {[—(5,0) U (^max, ^max + 5]} X [0, 

and the continuity conditions at the planar interfaces are 

E{0+,p) = E{0-,p), E,{0+,p) = E,{0-,p), 
E{Z^sj^+, p) = E{Z^a_^—, p), Ez{Z^i^^+^p) = Ez{Z,^s^^—,p). 



(19c) 



We shall sometimes find it convenient to adopt a general notation for both cases, by denoting 
the scalar transverse coordinate as = |a;^| and its domain by Q±. In the Cartesian case we have 
x± = x and i}± = [— ^max; -^max]; while in the cylindrically symmetric case we have x± = p and 
ni_ = [0 5 Pmax]- We shall also find it convenient to decompose the Laplacian as A — S^^ + A^, where 
A^ = dxx in the Cartesian case and A^ = ^dp{pdp) = <9p + j^dp in the cylindrically symmetric 
case. Physically, the transverse Laplacian term A±E leads to diffraction. 

Using this notation, the Cartesian system (jlSp and the cylindrically symmetric case system (jl9p 
are universally represented as 



E,,{z,x^) + A^E + k'o{u' + e\E\'nE = 0, 

^ ' (20a) 

{z,X^) G (0, Zmax) X 
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X±) G {[-(5, 0) U (Z,nax, ^max + 5]} X Ox, 

E{0+,x^) = E{0-,x^), Ez{0+,x^) = Ez{0-,x^), ^^^^^ 

^(^max+; = E(Zijii)_x — ,X±), Ez{Z-[n^^+, X±) = Ez{Zj^^:^ — , X±) . 

3.1.1. Local Transverse Boundary Conditions 

Following the approach first used in [i^], we set locally one-dimensional radiation boundary 
conditions of the Sommerfeld type in the transverse direction x^. To do so, we assume that the 
beam is localized around x± = 0, so that far from the beam center the nonlinearity becomes 
negligible, i.e., 

e\E\'^" < u"^, \x\ > Xmax or p> Pmax- 

Therefore, the field (approximately) satisfies the constant coefhcient equation: 

A£' + vlklE = 0, \x\ > Xmax or p> p^ax- 

We further assume that for \x\ > X^ax {p ^ Pmax) the field is composed predominantly of the 
outgoing plane (cylindrical) waves with nearly normal incidence on the boundary |x| = Xmarx 
{p = Pmax)- This leads to the following radiation boundary conditions in the 2D Cartesian case [201 ]: 

E^ - ikoi'oEl =0, Ez + ikoiyoE\ =0. (21a) 

In the 3D cylindrically symmetric case the local radiation boundary condition at p = Pmax reads (2l| : 

4-H'^\l'0kQPraa.yi) 

Ep-aE\^ =0, a = ^^, , 21b 

Hi) '(l^oKoPmax) 

where H^^ is the Hankel function of the first kind. The symmetry condition at the axis p = is 

^E{z,0)=0. (21c) 

We emphasize that these transverse boundary conditions are valid as long as the beam is localized 
around the axis and remains "far" from the transverse boundary at x = X^ax or p = Pmax- 

3.1.2. Nonlocal Longitudinal Boundary Conditions 

Similarly to the one-dimensional case (see Section [2. ip . the boundary conditions in the longitu- 
dinal direction z will be set in the linear regions at z = —6 and z = ^'max + They should render 
the boundaries transparent for all the outgoing waves, i.e., eliminate any non-physical reflections, 
and at the same time correctly prescribe the given incoming wave(s), see Figured! Unlike in the 
one-dimensional case, however, a two-way Sommerfeld boundary condition of type ([8]), which is 
local in the conflguration space, cannot be transparent for all the outgoing waves, because these 
waves travel with different longitudinal velocities that depend on their angle of incidence. 

Therefore, to accommodate all angles of incidence, we first separate the variables in the lin- 
ear Helmholtz equation (|20bp by expanding its solution with respect to the eigenfunctions of the 
transverse Laplacian. These eigenfunctions solve the ordinary differential equation: 

AxV^^)(xx) = -(A;J))VW, (22) 
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subject to the transverse boundary conditions (j2T|) . The resulting eigenvalue problem is not of 
the classical Sturm-Liouville type, since its operator is not self-adjoint (because of the radiation 
boundary conditions). As a result, the eigenfunctions are not orthogonal. Nevertheless, these 
eigenfunctions are bi-orthogonal [29|, Volume I] or, alternatively, real orthogonal, and still form 
a complete system. A comprehensive discussion on completeness of eigensystems arising in the 



diffraction theory, and on convergence of the corresponding series, can be found in 30(]. 

Since the system of eigenfunctions {-0^')} is complete, we can expand the field E and the 

fi z 
incoming beams E^^^ and -E^in"'"' as 



E{z,x^) = J2Mz)i^^^Hx±) 

(23) 

CXD OO ^ ' 



1=0 1=0 

In the transformed space, the linear Helmholtz equation ()20bp reduces to a system of uncoupled 
one-dimensional linear Helmholtz equations (ODEs): 

(^ + (fcf^)')n,(z) = 0, ikll^f = kl-{4^)', / = 0,l,...,oo. (24) 

Each of the uncoupled equations (j24|) formally coincides with equation ()6bp and has the same 
general solution composed of two waves one of which can be interpreted as propagation in the 
positive z direction and the other one — in the negative z direction. Unlike in equation ()6bp . 
however, the quantity {k^l^)'^ in equation (j24p may have a negative real part, in which case the 
waves become evanescent. It may also have a non-trivial imaginary part, which is due to the non- 
self-adjoint transverse (radiation) boundary conditions (see [20| for more detail). Regardless of the 
particular shape that the waves may assume, the longitudinal boundary conditions have to ensure 
that the field in the region z < —5 be of the form [cf. formula ([7])]: 

.,m 

u{z) = Uincie " +Cie II . 
Therefore, the two-way ABC at z = —6 can be written as 



dz 

Similarly, at the opposite boundary, z = Z^g^^ + S, we obtain: 



2ik\;^er<'yt:,l- (25a) 



Tz - • 



-2ikl!^e~''"\\^\^^^J. (25b) 



Boundary conditions ([25]) are local in the transformed space {ui{z)}^q. In this space, the two-way 
one-dimensional Sommerfeld conditions are applied independently for each individual mode defined 
by (|24p . The equivalent of relations (j25p after the inverse transformation of (1231) will result in a 
nonlocal pseudodifferential operator in the original space {E{z, x)}, see jiot ] or |3l| for more details. 
Therefore, the resulting boundary conditions are nonlocal two-way artificial BCs. 
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3.2. Discrete Approximation 

We build a semi-compact scheme for the Cartesian problem (jlSp . (j21ap . (j25p in Section [3.2.11 
and for the cylindrically symmetric problem (I19p , (I21b|) , ()21cP , (|25p in Section I3.2.2[ As in the one- 
dimensional case, we discretize the governing equations inside and outside the Kerr material, and 
then obtain a discretization at the material interfaces. The discrete transverse boundary conditions 
and the discrete two-way ABCs for both problems are described in Section 13.31 and Section 13.41 
respectively. 

3.2.1. 2D Cartesian Case 

On the rectangle [—3hz, Z^aa.^ + S/i^] x [— Xmax, ^max]) we introduce a uniform Cartesian grid of 
(iV + 7) X M nodes as 



Zn = n-h,, h, = ^, n = -3,-2,...,N + 2,N + 3, 

2X 

Xm = -Xma.^+ {m+l/2)ha:, hx = ^'"'^ , 771 = 0, 1, . . . , M - 1, 



(26) 



so that 

Zq = 0, = .^max) ^-1/2 — "^maxi ^M-l/2 — ^max- 

In this paper, we keep hz ~ so that all O (jiihx~''^ terms can be treated as terms of the same 

order k and denoted by O (h^\ For convenience, we also introduce the following notations for the 
field and the Kerr nonlinearity at the grid nodes: 

p 45f p/-^ ^ \ p 4^^ I p |2o-p 

Finally, we use the previous notation D for central difference operators, with the order of accuracy 
in the superscript and the differentiation variables in the subscript. For example, 

= - '^^r,^ + En,m-1 ^ ^^^^^^^ ^ ^ ^^,y 

Other notations for central differences are listed in Appendix [Bl 

To build a semi-compact approximation of the NLH (|18ap at the interior points n = l,...,A^ — 1, 
we first introduce the following mixed order discrete Laplacian: 

n(2) J? I n(4) p _ ^n~l,m. ~ 2En._m + -^'n+l.m 

I —En,m-2 + ISi^n.m-l — 30En,m + 16£'n,m+l ~ -£'n,m+2 ^^^^ 



12hl 

hi 4 

— ^-^n,m ~l~ ~^dzzzzEn.m ~\~ O {h ) . 

In order to remove the O {h?^ term on the right-hand side of (j27p . we consider the following 
expression that contains fourth-order derivatives with respect to both z and x, and approximate it 
to second-order accuracy using central differences: 

{dzzzz - dxxxx)En,m = ipzz " 9, J A£;„,,„ = [d^} - D^^]) AEn,m + O {h?) . (28) 
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Then, we employ the NLH (jl8ap itself and substitute the expression 



into formula (I28p . Next, we approximate the derivative dxxxxE in formula (I28p to second-order 
accuracy using central differences, and altogether obtain: 

dzzzzEn,m = -k^ (of) - D^^j) {u^Er,,m + ePn,m) + 4xL^n,m + O (h^) . (29) 



Substitution of (|29p into (j27p yields a semi-compact fourth-order discretization of the Laplacian, 
which leads to the following fourth-order scheme for the NLH (|18ap : 



^ zz ^ ^xx ^n.m ^xxxx j ^n,m 

+kl (l + ^Df) - (.2i?„,^ + eP.,^) = 0, 

n = 1, . . . ,iV - 1, m = 0,...,M-l. 



(30) 



To obtain a similar fourth-order scheme for the linear Helmholtz equation ()18bp . we repeat the 
previous derivation with ePn^m = and v = 1, which yields: 



. , n(2) I ( n{4) _ MM n(2) _ M n(2) \ , iu2 

Y2 I \^xx Y2^^^^ I 



-"n^m — >J) 



(31) 

n = -3, . . . , -1, -Fl, . . . , 3, m = 0, . . . , M - 1. 

Next, we consider material interfaces at the nodes n = and n = N. Using Taylor's expansion, we 
can write: 

— SSE'o/m + 108ii^i,m — 27£'2,m + 4£'3,.m 3/1^ „ „ „ /, 4\ 
— — — Oz-Ho+^m + -rr'^zz-'^O+.m + U [n 1. 

Don^ 11 ' 

Then, approximating the derivative dzzEo+,m with fourth-order accuracy: 

dzzEo+^m = A£^o+,m — dxxEo+,m = — ^o(^0+-^0,m + eO+-Po,m) — D^^J Eo^rn + C (^^) j 

we obtain: 

-85EQ rn + 108i?i m — 27£'2,m -|- 4i?3 ^ 



Deriving a similar formula for cJ^E'o-^m and equating the resulting expressions for (?^£'o+,m and 
dzEo-^rn, we get a fourth-order accurate approximation of the continuity condition Ez{0—) = 
Ez{0+): 

4^E-3^rn — 27£'_2,m + WSE^i^rri — 170£'o,m + 108£'i,m — 27i?2,m + 4£'3,m 



66hz 

H I 2 2 '™ y ^-'^xx -^o.m - U. (32j 
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Finally, substituting I'o-.m = eo-,m = 0, 1'o+.m = ^ and eo+,m = £) we have: 

'^E^S m — 27£'_2,m. + WSE^i m — 170£'o,m + lOSE'i m — 27S2,m + 4ii^3,m 



66/i 



z 



)6/z 
+ ^4t^Eo,m = 0. (33a) 



11 V 2 

A similar equation is obtained for the interface at n = A^: 

4-E'iV-3,m — 27£'7V_2,m + 108£'Ar_i_m — 170£'Ar,m + lOSE'Ar+i^m — 27EN+2,m. + 4i?7V+3,; 



66/i 

6h.k?: /I + _ 

4 
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^iV,„^ + ^PN,m^ + ^Di^jEM,m = 0. (33b) 



3.2.2. Cylindrically Symmetric Case 

We use the same grid (j26p . except that in the transverse direction we now have: 



{m + l/2)hp, hp = ^, m = 0,...,M-l, (34) 



so that 

P-l/2 = and Pm-1/2 = /'max- 

We also keep ~ hp so that all O (jiihp~-^^ terms appear of the same order O {h'^). 

To approximate the NLH (jl9ap at the interior points n = l,...,A^ — 1, we begin by introducing 
a mixed order discretization of the cylindrical Laplacian A = dzz + = dzz + "^p + ^i9p : 



(Og^ + + —D^;A En,m = ^En^m + ^dzzzzEn,m + O (h^) . 

\ Pm / -Lz 



(35) 



To remove the O (/i^) term on the right-hand side of (|35p . we start with the second-order 
central difference approximation of the expression {dzzzz — ^'ji)En,m = {dzz — ^p)^En,m, where 
Ap = p~^dp — p~^dpp + 2p~^dppp + dpppp , and using the NLH (|19ap itself, obtain: 

dzzzzEn,m = - ( Of) - Of) - —Df~^\ {l^^E^.m + ePn,m) 

\ Pm J (36) 

+ (o~^ - + 2o^^ , r)(2) \ p ,0 /, 2^ 

' \Pm Pm ^pp ^ ^Pm ^ ppp ^ ^ pppp I ^n,m -T ^ \n ). 



Substitution of (j36j) into (j35|) yields a semi-compact fourth-order discretization of the cylindrical 
Laplacian, which leads to the following fourth-order scheme for the NLH ([T9a]): 



( I n(4) , J_ n(4)^ E (o-^D'^'^^ - o-^D^^) , On"! n(2) I n(2) ^ f 

\^zz ^ ^pp ^ p ^ p J ^n,m ^2 \Pm ^ p Pm ^ pp ^ ^Pm ^ ppp ^ pppp J ^n,m 



V Pm / _ 

n = 1, . . . , - 1, m = 0,...,M-l. 
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To obtain a similar fourth-order scheme for the hnear Helmholtz equation ()19bp . we repeat the 
previous derivation with e = and v = 1, which yields: 



( L>(2) I 2^(4) , J_I)(4)^ E (o''-^ _ o-^^m I 2n-l Z)(2) + Z){2) E 

\^zz ' ^pp ^ p ^p J '^n,m ^2 \ym ^p Pm. ^ pp ^ ■^Hm ^ ppp ^ PPPP J ^ri,m 



1 + Mfl)(2)_^(2)_±^(2) 



n = -3, . . . , -1, iV + 1, . . . , 3, m = 0, . . . , M - 1. 



F - n 



The analysis of material interfaces at n = and ?i = is very similar to that of Section 13.2.11 
and we arrive at the following fourth-order accurate approximation of the continuity condition 
£;,(0-) =E,(0+): 

4-E^-3,m — 27£'_2,m, + 108£'_i^m, — 170£'o,m + lOBE'i^m — 27£^2,m + ^E^,m 

Q6hz 

( + 5z±i£p„, J + ^ ( o'S + ^D'*i] = 0. (39) 



11 V 2 2 
Substituting i^o-,m = 1, eo-,m = 0, fo+,m = and eo+,m = e into ([39|) . we have: 

4-E-3,m — 27-E'_2,m + 108i?_i^m " 170£'o,m + lOSE'i^m " 27£^2,m + ^E^,m 

Q6hz 

A similar equation is obtained for the interface at n = A^: 

4-E'Af-3,m - 27£'jV„2,m + 108-E'jV-l,m — 170-5' jV,m + lOS-EAr+l.-m " 27£^Af+2,m + 4£'jV+3,m 

+^ (^i?^,,n + + '-^ {d(;J + E^,^ = 0. (40b) 

5.5. Local Transverse Boundary Conditions 

In this section, we briefly describe a discrete approximation of the transverse boundary con- 



ditions (j2ip . In doing so, we follow the approach of j2l|, where additional details can be found. 
Let us first consider the radiation boundary conditions ()21ap and (|21b[) at the "upper" boundary 
m = M — 1/2. We will use their discrete counterparts to express the values of the field at the 
ghost nodes En^M and E^^m+i via the values at the inner nodes En,M-3, and En,M-i, and 

thus eliminate the ghost nodes. A fourth-order approximation of either Cartesian or cylindrical 
radiation boundary condition centered around m = M — 1/2 (which corresponds to x = Xmax or 
P = Pmax) is given by 

En,M-2 — '27En,M~l + 27En,M — En,M+l 

24h± 

— En,M-2 + ^En,AI~l + 9En,AI — -E„ A/+1 

re = °' 
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where in the Cartesian case a = i^oko, see formula (j21ap . and in the cyhndrical case a is defined 
in (|21bj) . see 



211]. Equivalently, we can write: 

[c_2, 



■n,M-2 



E; 



'n,M+l 



0, 



where 



[c-2, 



[1, -27, 27, -1] 



2ah 



[-1, 9, 9, -1] 



However, specifying this boundary condition alone is not sufficient, because the fourth-order finite 
difference equation that we use in the x± direction requires an additional boundary condition. The 
choice of the latter allows for more flexibility as long as the resulting method is fourth-order accurate 
and stablelfl Hereafter, we choose this second condition as the fourth-order accurate extrapolation 
of the ghost value E^^m+i via {En,M-2,^ ■ ■ ■ ■, En,M}, which can be conveniently written as 



E, 



n,M+l 



'n,M+j- 



Combining the two discrete boundary conditions as 



0, C_2, C_i, Co, Cl 

-14 -6 4 -1 



E, 



n,A/-3 



_En,M+l_ 

we can express the ghost values En m and En m+i in terms of the interior values: 



En,M 
En,M+l 



1 



Co + 4ci 



— Cl c_2 + 4ci c_i — 6ci 
Co 4c_2 - 4co 4c_i + 6co 



En,M-3 
En,M-2 
En,M~l 



(41a) 



In the Cartesian case, the derivation is repeated to obtain the discrete discrete radiation boundary 
condition at x = — ^max, i-e., at m = 0: 



En 

E.. 



n,-l 



1 



Co + 4ci 



4c„i + 6co 
c_i — 6ci 



4c_2 - 4co 
C_2 + 4ci 



Co 

-Cl 



Enfi 
En, I 
En,2 



In the cylindrical case, the symmetry (j21cp is enforced as follows: 



E. 



n,-l 



E, 



n,0, 



E, 



n,l- 



(41b) 



(42) 



Note also that in the Cartesian case there is an alternative way of building the discrete transverse 
boundary conditions. It does not require a finite difference approximation of the continuous bound- 
ary conditions (j21ap . and is rather based on analyzing the roots of the fourth-order characteristic 
equation that corresponds to the five node discretization in the x direction. The idea is similar to 
that behind boundary conditions (|17p . and the reader is referred to for more detail. 



Stability of these approximations can be studied by the methodology of [32 
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3.4- Nonlocal Longitudinal Boundary Conditions 

In this section, we construct a discrete counterpart for the two-way ABCs (j25p . In the continuous 
case of Section 13.11 we separated the variables in the hnear Helmholtz equation (|20bp outside the 
Kerr region, and then obtained the ABCs in the transformed space. In the discrete case, we also 
begin by separating the variables in the Cartesian (|3ip and cylindrical (|38p difference Helmholtz 
equation at the exterior grid nodes: 



m 



0,, 



M-h 



n 



0,-1, 



n 



N,N + 1,N + 2. 



Subsequently, we derive the ABCs in the transformed space. This derivation is identical for the 
Cartesian geometry of Section 13.2.11 and the cylindrical geometry of Section 13.2.21 

We first identify the transverse components in the finite difference operators of (j3ip and ()38p . 
The transverse part of the discrete Laplacian for the Cartesian case is 



£,(4) 



12 "^"^ 



Mn(2) 

-j^2 i 



(43a) 



whereas for the cylindrically symmetric case it is given by 



klhl 
12 



pp 



rm ^ o 
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PP n P 



(43b) 



-2^(2) , 2n-lD{2) + I)(2) 
rm pp ' ^rm ^ ppp ' ^ pppp 



The separation of variables in equations (j3ip and (|38p will be rendered by expanding the solution 

T 



with respect to the transverse eigenvectors "0^'^ 
satisfies the following difference equation on the grid [cf. equation ([22 



LV(^) = -(fc«)V«, 



0,1,, 



M- 1. 



Each eigenvector -i/^^'-* 



(44) 



In the Cartesian case, the operator in (I44p is defined by formula (I43ap . and the solution '^^^ 
is subject to boundary conditions (|41ap . (|41bp . In the cylindrically symmetric case, the opera- 
tor L"*- in (j44p is defined by formula (j43bp . and the solution -f/'^^) is subject to boundary condi- 
tions (j41ap . (j42p . The argument behind linear independence of l^*^')} in the Cartesian case is 
based on bi-orthogonality (real orthogonality) of the eigenvectors and can be found in [i^]. For the 
cylindrically symmetric case, the continuous eigenfunctions are also real orthogonal, but the discrete 
eigenvectors are not, see [2l|. Yet we observe numerically that they are linearly independent. 
The M linearly independent eigenvectors are convenient to arrange as a column matrix: 



dcf 



(0) 



(M-l)- 



that will diagonalize the discrete transverse Laplacian, i.e., L^"^ = ^'A, where 

A=dlag{-(fcf)^-(fcW)^...,-(fcf-^))2}, 
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and the eigenvalues — {k±^^ ^^'^ defined in 

It will also be convenient to consider the following M-dimensional vectors: 



that contain the values of the field arranged in the transverse direction. With this notation, we can 
recast both scheme (1311) and scheme (1381) in the vector form: 



hi 

For each n, let us introduce the vector variable 

U = ^~^E 

^ ri ^ *-'r). 



0. 



(45) 



(46a) 



Equality = is the expansion of £n with respect to the eigenvectors ip^^^ , where the coefficients 
are given by the components of Un- Similarly, we can expand the incoming beam profiles: 



TjO =^-ir.o 



^inc ^ ^inc J 



u: 



4£f 



u 



inc,0 



U. 







Zmax dcf 



inc,M-l. 

Zmax 
^inc,0 



Zmax 
inc,M-l. 



,0 



Zmax dcf 



-^inc,0 



.-^inc,M-l. 

■ piZmax 
-'^inc.O 



iTiZmax 
.-^inc,Af-l. 



(46b) 



(46c) 



Formulae (j46p are discrete counterparts of (j23p . Substituting expansions (j46p into equation ([45 
and diagonalizing L"*-: L^Sn = L-^^Un = ^'ACZ„, we have: 



(47) 



, Un M~i]^ , we obtain: 



Finally, multiplying equation (j47p by the inverse matrix ^ ^ from the left we separate the variables. 
Recasting the result via individual components of C/„ = [un,o,Un^i 



1 + 



12 



hi 



(48) 



/ = 0,1,...,M-1. 



Formula (|48p is a system of M uncoupled ordinary difference equations, which is a discrete coun- 
terpart of the continuous uncoupled system ()24p . 

Each of the uncoupled difference equations (j48p is identical to the one-dimensional difference 
equation (fTSj) if we redefine A;^ of (fT5]) as 



(0^2 



k 



2 _ /cq (A;_|_ ) 



{k 



(0x2 
II ^ 



i-2/,2 
12 



12 
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Therefore, similarly to (jl7ap we can write for the ghost node n = —4: 



U-A,i = {qi - qi) Qi u^n^ i + qiu_sj 



where qi and qf^ denote roots of the characteristic equation for a given I, and the incoming compo- 
nents u?^^ I are defined in (j46p . Recasting the previous equality in the matrix form and transforming 
back into the configuration space, £ = "ifU, we obtain the two-way discrete ABCs: 



Qo 



^-^£l, + ^ 



qM~i 



(49a) 



Likewise, for the ghost node n = -|- 4 we write similarly to (jl7b|) : 

UN+4,l = {qT^ - qi) ^r^^^cj + QlUN+3,1, 

and arrive at the following two-way discrete ABCs: 



£. 



N+4 



gp ^ - go 

ql 



go 



gA/-l 



N+3- 



(49b) 



Relations (j49p provide a fourth-order accurate approximation to the boundary conditions (j25p for 
5 = 3hz- In the Cartesian or cylindrical case, relation (j49ap is substituted into equation (j3ip or (j38p . 
respectively, for n = —3, and relation (I49bp is substituted into equation (f3T]) or ([38l) . respectively, 
for n = A -|- 3. This eliminates the ghost values from the schemes (j3ip and (j38p and thus closes 
the system of finite-difference equations on the grids (j26p and (|34p . 



5.5. Extension to the Multi-Layer Material 

In the case of a grated Kerr material described in Section [1.31 there are additional discontinuity 
points defined by formula (j5ap . The interface conditions at each discontinuity point z are the same 



as at z = and z = Zrr 



E{z+,x±) = E{z-,x±), 



dE , 
dz 



[z+,x±) 



dE 

dz 



(z-,x_l). 



In the simple case when z coincides with one of the grid nodes, the discrete approximation of the 
continuity conditions is given by the same formula as (j32p for the Cartesian case and by the same 
formula as (j39p for the cylindrical case. Hence, to solve the one-dimensional NLH for the multi-layer 
case one needs to apply the corresponding discrete interface condition of type ([32|) or ([39]) at each 
plane (j5ap . The extension to the case when a discontinuity plane does not coincide with any of the 
uniform grid surfaces ([26]) or ([M]) can be obtained by building separate grids for different layers. 
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4. Newton's Solver 



Here we briefly outline our approach to building a Newton type solver for the NLH. The reader is 
referred to [l^ . Section 3] for a detailed description. The schemes for the NLH that we constructed 
in Sections [2] and [3] lead to systems of nonlinear difference equations that we symbolically write as 

F{E) = 0, 

where the quantities F and E are complex. In the one-dimensional case, they are vectors of 
dimension + 7, which is the dimension of grid 1^, F, E G C^^''. In the two-dimensional case, F 
and E can be interpreted as matrices of dimension + x M, which is the dimension of grid (I26p . 
and for convenience we reshape them as (A^ + 7)M-dimensional vectors: F, E & f£{N+7)M ^ 

To linearize the transformation F{E), we first notice that the Kerr nonlinearity P = \E\'^'^E 
is Frechet nondifferentiable as long as E is complex. To overcome this, we separate the real and 
imaginary parts and recast the field E and mapping F as real vectors of twice the dimension: 

F{ E) : C(^+^)^^ ^ c(^+7)M _^ p^^y_ j^2(7V+7)M ^ ^2iN+7)M ^ 

The new transformation F{E) is differentiable in the conventional real sense. Differentiation results 
in Newton's linearization of F that involves the Jacobian J , and leads to Newton's iterations: 



1 

dcf 



J{E 



0> 



F{E^^'^). (50) 



The convergence of Newton's method is known to be very sensitive to the choice of the initial 
guess. In our experiments, we take the simplest initial guess E^^^ = 0. We have also observed 
numerically that the algorithm was more likely to converge if, during the first stage of the iter- 
ation process, when the iterations E^^^ are "far" from the solution, we introduce the relaxation 
mechanism: 



max{l,||5i;0+i)||oo} ^ ' 

where w G (0,1]; typically uj = 0.5. While this mechanism enables the algorithm to converge 
for a wider range of cases, it also slows down the convergence (from quadratic to linear rate). 
Therefore, once the iterates E'^^^ are "sufficiently close" to the solution so that ||oo < 0.01, 

we change back to w = 1, thereby reverting to the original Newton's method (|50|) . The criterion 
for convergence that we employ is the inter-iteration distance threshold jJ^/'^-'-'l < 10~^^. 



5. Summary of the Numerical Method 

The NLH ()20p subject to local transverse boundary conditions (j2ip and nonlocal two-way lon- 
gitudinal boundary conditions ([25]) is discretized on the grid (126]) or (fM]) . In the Cartesian case, 
we obtain semi-compact schemes (j30p and (|3ip in the interior and exterior of the Kerr material, 
respectively, and discretization ()33p for the continuity conditions at the interface. In the cylin- 
drically symmetric case, we arrive at the semi-compact schemes dSTj) and (|38p at the interior and 
exterior nodes, respectively, and discretization ()40p for the continuity conditions at the interface. 
For both geometries, we also employ discretization (|4ip for the local transverse radiation boundary 
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condition, and discretization (j49p for the non-local two-way boundary conditions at z = —3hz and 
z = Zj^i^^ + 3hz- In addition, discretization ()42p is used at the axis of the cylindrical system. 

The resulting system of nonlinear difference equations with respect to complex unknowns En^rn 
is recast in the real form at the expense of doubling its dimension. Then, Newton's linearization 
is applied, see formula (jSOp . which yields a 2{N + 7)M x 2{N + 7)M sparse Jacobian matrix, with 
the bandwidth of 2M for the interior and exterior grid points, where n 7^ 0, A^. For the points at 
the interfaces, where n = or n = N, the bandwidth is 6M. 

At each Newton's iteration (j50p or (jSip . this Jacobian needs to be inverted. Currently, we are 
using a sparse direct solver to invert the Jacobians. This entails an O • M^) memory cost and 
hence imposes a fairly strict limit on the grid dimension. For example, a typical grid dimension of 
N X M = 1000 X 320 results in the memory requirement of about 6Gb. 



6. Finding an NLS-Compatible Incoming Beam 

As indicated in Section 11.11 the NLH is the simplest nonparaxial model that generalizes the 
NLS. Accordingly, one of our key goals is to investigate how the addition of nonparaxiality affects 
the solution. In order to do so, we shall use the NLS solutions as "benchmarks," and compare them 
with NLH solutions computed for "similar" input parameters. 

We note that for an incoming beam E^^^^^ which impinges on the material interface at z = 
0— , only a part of it that we denote by -Brcfractcd passes through whereas the rest gets reflected. 
In contradistinction to that, in the NLS framework all of the incoming beam -E^J^^^^ propagates 
forward. Therefore, in order to have comparable incoming beams for these two models, we should 
choose the NLH incoming beam E^^^^ so that the refracted part of it at 2; = 0-|- be close to the 
NLS initial data -EJ^^^^, i.e., 

sSted(0+, x^) « E^'^'ix^). (52) 

A comprehensive solution to this problem is nontrivial, because the reflection at the nonlinear 
interface z = depends on the NLH solution itself for z > 0. Therefore, in this paper we use an 
approximate treatment which experimentally proves sufficient. 

In order to present this approximate treatment, let us first consider the one-dimensional linear 
problem: 

—\l + ,y\z)E = 0, i^{z) = { ' (53a) 

dz^ z > 0, 

with the wave E^^^^e^^ impinging on the interface from the left. The overall field has the form: 

lET''e- + Re-, -oo<.<0, 
[Te'"^, < 2 < 00, 

where R and T are the reflection and transmission (refraction) coefficients. The values of R and T 
are obtained from the continuity condition at the interface 

E{0-) = E{0+), ^(0-) = ^(0+), (53c) 



which yields: 



l«cted(0+)| = |T| = ^|i^r"^|. (54) 
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Formula (j54p is a standard result for the transmission of plane waves with normal incidence at a 
single linear interface, see, e.g., [s^, Section 7.3, eq. (7.42)]. 

We shall use this simple refraction formula to approximate the refracted beam of our weakly 
nonlinear multi-dimensional problem: 

£^Sted(0+,^^) - , ' 

Next, we assume that the NLH solution is close to the refracted incoming beam, E{Q+,x±) ~ 
-E^refracted(0+, x^), and obtain: 

^Sted(0+,a.±) - 1 ' ^^""{^l-)- 

1 + + e |Srefractcd(0+, X^)\^'' 



Finally, requirement ()52p implies: 



pO,NLS/ 



Equation ([55]) will be used throughout Section [7] for all collimated incoming beams. 



7. Numerical Experiments 

7.1. 2D Cubic NLH (Solitons) 

7.1.1. A Single Collimated Beam (Nonparaxial Soliton) 

The Cartesian configuration (D = 2) models propagation in planar waveguides. In the case of 
a cubic nonlinearity (a = 1) and = 1, the one-dimensional NLS (j4]) has solitary wave solutions: 

E(z x) = ^'^ exp {ik^zjl + 4)) ^ V2 exp(ifcoz(l + i(fcoro)-2)) ^^^^ 

' V ^ / cosh(/A;oa;) kor^^ cosh(2;/ro) 

which are called solitons. In formula (j56p . tq is the soliton width and / = {k^ro)^^ is the nonparax- 
iality parameter, which can also be interpreted as the reciprocal beam width measured in linear 
wavelengths: 27r/ = Ao/tq. 

We solve the Cartesian NLH (|18p on an elongated domain: .^max = 240, Xmax = 12, and for 
/co = 27r/Ao = 4, = 1, and e = k^"^ . The problem is driven by the incoming beam 



l-hA/l + esech^ (x/V2) / 
El^{x) = V _ g^^j^ (■^/^ 

for which the refracted beam is (approximately) an NLS soliton profile of the width tq = ^/2: 

-^refracted ^ ^ech (x/a/2) , 



see formula ([55]) . The corresponding nonparaxiality parameter is / = l/\/32 w 0.177, which means 
ro = 0.90 • Ao and which is considered a very narrow beam. 
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Figure 3: (color online) A nonparaxial soliton for the 2D Cartesian NLH with cubic Kerr nonlinearity. A) Normalized 
on-axis |_Ep (black, dotted) and Sz (red). B) |_E|^ contour plot (zoom in on several oscillations). C) Sz surface plot. 

In this simulation, the field was assumed symmetric with respect to the x-axis, E{x) = E{—x). 
This allows us to increase the resolution in the x direction by a factor of two. A non-symmetric 
simulation at half the resolution provides very similar results. The grid dimension that we took was 
N X M = 4480 X 112, which translates into the resolution of Xo/hz = 30 and Xo/hx = 15, i.e., 30 
grid points per linear wavelength in the z direction (axial) and 15 grid points per linear wavelength 
in the x direction (transverse). 

In Figure [3]A., we plot the on-axis amplitude of the Cartesian NLH solution. The square ampli- 
tude \E\'^ exhibits fast oscillations in the z direction, as can be seen in the insert of Figure[3]A, and 
in Figure [3jB. Although at a first glance these oscillations may appear a manifestation of numerical 
instability, in fact they are physical and indicate the presence of a backward propagating component 
of the field. Indeed, for a field with both forward and backward propagating components: 



(57) 



the square amplitude is given by the expression: 

|Sp ^ \A\^ + 2 Re (^AB*e^''"'^^ + (58) 

which has a ~ 2ko oscillating term. We note that the amplitude oscillations in Figures [3]A and [3(3 
are indeed ~ 2kQ, as predicted by formula ()58p . We further note that these oscillations are also 
exhibited by the explicit solutions of the ID NLH 0]. 

We recall that for the NLS the square amplitude is proportional to the energy fiux densitylf] 
and that the L2 norm of the solution H^Hl = / is a conserved quantity proportional to the 

total energy fiux or, equivalently, the beam power. For the NLH, however, the proper measure of 
the energy flux density is the Poynting vector: 

S = fco^Im {E*VE) 



^ In the Gaussian system, the quantity has the units of energy flux: 

per unit of time. 



cm^ - sec 



i.e., of energy per unit area 
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rather than the square amphtude. Accordingly, the conserved beam power (i.e., the total energy 
flux) is the integral its z component over the beam cross-section: 



N 



Then, for the field ()57p with both forward and backward propagating components, the energy flux 
density reduces to 

S,^{\A\'-\B\^), (59) 

i.e., to the difference of the forward and backward square amplitudes. Clearly, expression (j59p 
contains no (rapidly) oscillating terms. The Poynting flux Sz for the 2D NLH solution is given 
in Figure [3lC, and is indeed much smoother than the amplitude, see Figure [3K. We therefore 
suggest that the energy flux density provides a more adequate quantitative measure of the long- 
scale (collapse) dynamics in the NLH model. 

The key physical question that the simulations in this section attempt to answer was whether 
there exist any solitons beyond the paraxial limit, i.e., of the 0(Ao) radius. Considering the energy 
flux of the 2D NLH solution with a = 1 shown in Figure [3p, we see that it indeed resembles a 
soliton propagating essentially unchanged in the positive z direction. We can therefore conclude 
that such a nonparaxial soliton does exist. 

Let us also note that nonparaxial solitons (solutions of the NLH, rather than NLS) for a single 



collimated beam were obtained in [3J] for the case of a semi-inflnite Kerr medium. Our formulation 



is different as it involves a flnite-width Kerr material slab with the interfaces that may partially 



reflect the waves. Hence, a direct comparison of our results with those of [3j] is not appropriate. 



However, a comparison from the standpoint of physics may be of interest for the future. 

7.1.2. Grid Convergence Study 

In order to demonstrate the fourth-order grid convergence in the nonlinear regime, we simul- 
taneously reflne the grid in the transverse and longitudinal direction, and monitor the maximum 
difference between the computed flelds for each pair of consecutive grids, the coarser and the finer, 
that differ by a factor of 2 in size. For the grids with fewer than roughly 7 points per linear 
wavelength, the iterations diverge, apparently due to insufficient resolution. Hence, we choose 
our coarsest grid to have Xo/hz = 7.5 nodes per wavelength, and compare the results with those 
on the twice as fine grid, Xq/^z = 15. Then, we keep decreasing the size and hence increasing 
the dimension of the grid, and the largest dimension that we can take is limited by the memory 
requirements of the LU solver that we employ for inverting the Jacobians (see Sections H] and [5]) . 
Currently, it is close to x M = 4480 x 112, which corresponds to 30 points per linear wavelength 
in the z direction. The results of the grid convergence study are summarized in Table [TJ The 
convergence rate that we find is O (/i^'^) , which is close to the design rate of O (/i^) . 



Table 1: Grid convergence study for the 2D Cartesian NLH with a = 1, ko = A, e = kg^ , Zmax = 240, and Jfmax ~ 40. 



{hz,hp) 


[Xo Ao\ 


(Xo Ao\ 


(Xo Xo\ 


Vl5'7.5y 


V2I' 10 J 


V30' 15 J 


11^(2/.) 


1.1 


0.30 


0.080 


log2||i?(2^^-i?W||oo 


0.20 


-1.7 


-3.6 
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7.1.3. Collision of Nonparaxial Solitons 

The NLH is an elliptic equation with no preferred direction of propagation. Therefore, it can 
be used to model the interaction of beams traveling at different angles, and specifically counter- 
propagating beams. To demonstrate this capability, we solve the 2D NLH with a = 1 for two 
configurations. In the first one, two perpendicular nonparaxial solitons collide, while in the second 
one, two counter-propagating beams collide almost head-on, at the angle of 150°. Note that the 
paraxial approximation is invalid in the region of interaction between the beams for either case. 
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Figure 4: (color online) Collisions of two nonparaxial solitons in the 2D Cartesian NLH with cubic Kerr nonlinearity. 
A) Collision angle 90°, Sz surface plot. B) Collision angle 90°, surface plot. C) Collision angle 150°, Sz surface 
plot. D) Collision angle 150°, \E\'^ surface plot. 



For the perpendicular beam configuration, we solve the 2D NLH with = 6, Z„ 



20, 



-^max 

at Z = 



= 30, 
0, X 



1, and e = k. 







The forward-traveling incoming beam enters the material slab 



10, and propagates in the — 45° direction, while its counterpart enters at z — 



X = 10, and propagates in the —135° direction. The resolutions were Xo/hz = Xo/hx = 10 points 
per linear wavelength. A surface plot of the energy flux density Sz is shown in Figure HJA.. Positive 
values of Sz (forward propagation) are red, while negative values (backward propagation) are blue. 
As in the paraxial NLS model, the two nonparaxial solitons are almost unchanged by the collision. 
A surface plot of \E\'^ is shown in Figure HB; the oscillations in the interaction region are due to 
the presence of counter-propagation waves. 

For the head-on collision configuration, we solve the 2D NLH with ko = 4, 2max = 30, X^ax = 
12, = 1 and e = The forward-traveling incoming beam enters the material slab at z = 0, 
X = 4, and propagates in the —15° direction, while its counterpart enters at z = .^max, x = 4, and 
propagates in the —165° direction, resulting in a collision at the angle of 150°. The resolutions 
were Xo/hz = Xo/h^ = 16 points per linear wavelength. The results presented in Figures [DC and 
HP show that as in the previous case, the solitons re-emerge essentially intact after the collision. 
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1.2. Arrest of Collapse in the NLH 
7.2.1. 3D Cylindrically Symmetric Case 

We solve the cylindrically symmetric NLH (fTOl) for a = I, ko = 2'k/\q = 8, u = 1, ^max = 9, 



/ 2 — ^ 

and pmax = 3.5. The problem is driven by the incoming beam Efj^^{p) = ^^"^ ^1^'^'^ — —e"'^ , for 
which the refracted beam is approximately a Gaussian: -E'^'cfj-actcd ~ 7 see formula ([5^ . The 
grid dimension is x M = 1080 x 360, which translates into the resolution of Xo/h^ = 94 and 
Xo/hp = 81, i.e., 94 grid points per linear wavelength in the z direction (axial) and 81 grid points 
per linear wavelength in the p direction (radial). 

While this estimate shows that the waves in the linear region are very well resolved, we note 
that the NLH 

AE + 4l {\E\^) E = 0, 4l = ^0 (1 + el^n , 

supports waves with nonlinear wavenumber /cnl- In order to ensure that these nonlinear waves are 
also well resolved, a similar resolution estimate should be performed for the nonlinear wavelength 
Anl = Ao/\/l + Below, we will see experimentally that at the maximum self focusing point 

(with the maximal amplitude), we have e\E\'^ ~ 4.6. Hence, the nonlinear waves with the minimum 
wavelength of Anl ~ Ao/2.4 are well resolved, with Anl/^z ~ 40 points per nonlinear wavelength 
in the z direction, and Anl/^p ~ 31 points per nonlinear wavelength in the p direction^ 

The nonlinearity coefficient was e = 0.15. The parameter that controls the beam collapse in the 
corresponding critical NLS <^ is the ratio of the incoming beam power Pq = pe~'^^ dp = ^ to 
the critical power Pc ^ 1.8623/(e/cg), see For the NLH with the values of the parameters 
we have chosen, this power ratio is related to the nonlinearity coefficient e as 

p= — Ki ko = 1.29. 

^ Pc 4-1.8623 

In Figure OA, we compare the cylindrically symmetric NLH solution with the corresponding 
NLS solution at the axis of symmetry p = 0. Since the beam power is 29% above Pc, the solution 
to the NLS blows up and its on-axis amplitude tends to infinity at ^; ~ 5.5. The corresponding 
NLH solution, however, remains bounded and its amplitude attains its maximum max\En m\ ~ 5.5 

n.m ' 

at z ~ 6.25. This yields the maximum Kerr nonlinearity of maxjelE'nmP} ~ 4.6. 

n,m. ' 

The square amplitude and energy flux density of the cylindrically symmetric NLH solution 
are displayed in Figure [6l As in the "soliton" case, fast oscillations in the z direction are clearly 
observed for the square amplitude, but not for the energy flux, which appears smooth. 

7.2.1.1. Crid Convergence Study. In order to demonstrate the fourth-order grid convergence for 
the cylindrical geometry case, we conduct a grid convergence study similar to that of Section [7.1.21 
For the grids with fewer than roughly 18 points per linear wavelength, the iterations diverge. This is 
apparently due to insufficient resolution in the region of strong focusing, where it will only be about 
18/2.4 = 7.5 nodes per wavelength. As such, the coarsest grid we have taken had 17.5 points per 
linear wavelength in the z direction, and the finest grid was N x M = 1140 x 380, which corresponds 
to 100 points per linear wavelength in the z direction. The results of the grid convergence study 



For the soliton simulations in Section Pm the nonlinearity was smaller and Anl « Ao , so that a separate resolution 
estimate for Anl was not needed. 



29 



40 




A B 



Figure 5: (color online) Arrest of collapse in the NLH: Normalized on-axis square-amplitude \E\^ (blue solid), the 
Poynting vector Sz (red dashed), and the NLS solution (black dotted) on the axis. A) D = 3 and a = 1. B) _D = 2 
and a — 2. 



are summarized in Table [2j The convergence rate that we find is O (/i^'®^), which is even somewhat 
better than the O (/i^) theoretical rate. 

Table 2: Grid convergence study for the cylindrically symmetric NLH with a — 1, p = 1.29, Zmax ~ 9, and pmax = 3.5. 



(hz, hp) 


/Ao Ao\ 
V35' 30 J 


/Ao Ao\ 
V50' 43 J 


/Ao Ao\ 
V70' 60 J 


/ Ao Ao\ 
VlOO' 86 J 




3.63 


0.965 


0.176 


0.0225 


log2||i?(^'^^-i^W||^ 


1.86 


-0.051 


-2.51 


-5.47 



1.2.1.2. Effect of the Domain Size. Our simulations show that the convergence of Newton's iter- 
ations depends on the domain size, and specifically on the length .^max of the Kerr material slab. 
To investigate this dependence, we attempt to solve the (2 + 1)D NLH for several domain sizes 
^max = 1; 2, 3, ... , 15. In order to limit possible effects of the transverse boundaries on the conver- 
gence, it was positioned relatively far from the axis, at pmax = 5. In order to limit possible effects of 
under-resolution, we have chosen moderate resolutions of Xo/h^ = 53 points per linear wavelength 
in the longitudinal z direction and Xo/hp = 31 points per linear wavelength in the transverse p 
direction. The results are displayed in Table [3l It can be seen that for some domain lengths the 
algorithm converges, while for others it diverges. It may be possible that the divergence observed 
for the domain lengths between 4 and 7 is related to the boundary z = .^max being positioned too 
close to the region of maximum self-focusing, see Figure [5]A.. 
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Figure 6: (color online) Arrest of collapse in the cylindrically symmetric NLH. Plots of the square amplitude (top) 
and the energy flux density (bottom). 



Table 3: Convergence of Newton's method for the cylindrically symmetric NLH with a = 1, fco = 8, and e = 0.15 
on the series of domains with Zmax = 1, 2, . . . , 15 and pmax ~ 5. The criterion of convergence is < 10^^^ and 

to = 0.5. 



^max 


1-3 


4-7 


8, 9 


10-15 


Convergence 


YES 


NO 


YES 


NO 



7.2.2. The 2D Quintic Nonlinearity Case 

We solve the Cartesian NLH (fT8]l for a = 2, ko = 2tt/\q = 8, = 1, 2max = 6, and Xmax = 3. 

/ A 2 

The problem is driven by the collimated incoming beam 

refracted beam is approximately a Gaussian: E^c^^.^x.cA ~ see formula ([25]) • The grid dimension 
is X M = 900 X 300, which translates into the resolution of A0//12 = 120 grid points per linear 
wavelength in the z direction and Xo/hx = 80, grid points per linear wavelength in the x direction. 
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Figure 7: (color online) Arrest of collapse for an inclined beam in the 2D Cartesian NLH. Plots of the square amplitude 
(top) and the energy flux density (bottom). 



The shortest nonhnear wavelength was Anl = Ao/\/l + emax \ E\^ ~ Ao/1.85. The nonhnear waves 
are therefore still well resolved, with X^i^/hz = 65 and Anl/^x = 43. The nonlinearity coefficient 
was chosen e = 0.125, and the ratio of the incoming beam power was Pq/Pq ~ 1.30. The results 
are displayed in Figure [SJ3, and are similar to the 3D cylindrically symmetric critical case. 



7.2.3. An Inclined Beam: Focusing- Defocusing Oscillations 

We solve the 2D NLH with a = 2, ko = 8 and z^^ = 1 on the domain with ^max = 12 
and Xinax = 12 for a Gaussian incoming beam entering the Kerr material at z = 0, x = 4 and 
propagating at the angle of -7r/5.3. The nonlinearity coefficient was e = 0.12, which yields the input 
power of 28% above critical. The grid was N x M = 400 x 800, which corresponds to resolutions 
of Xo/hz = Xo/hx = 26 points per linear wavelength, and 14 points per nonlinear wavelength Anl- 

As shown in Figure O the beam undergoes two focusing-defocusing oscillations, which qualita- 
tively agrees with the predictions of the modulation theory for the NLS [171]. This is the first time 
that two focusing-defocusing oscillations are observed in a critical NLH model. 
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1.3. The Effect of Adjusting the Incoming Beam 

As indicated in Section [6l the incoming beam for the NLH needs to be adjusted so that to 
enable a more accurate comparison of the results with those obtained for the corresponding NLS. 
In this section, we investigate the difference between the NLH solutions obtained with or without 
adjusting the incoming beam. Namely, we analyze the critical case D = 3, a = 1, and rerun the 
simulation of Section [7.2.11 with .^max = 8.5 and for the incoming beam E^^^ = e~P , i.e., without 
adjusting the incoming beam. The resolutions are Xq/Hz = 83 and Xo/hx = 67 points per linear 
wavelength. The results presented in Figure [8] show that in this case the collapse occurs later 
and achieves a smaller maximum self-focusing than for the adjusted incoming beam. The insert of 
Figure [8] also shows that near the boundary (after the refraction by the interface) the solution with 
the adjusted incoming beam is indeed much closer to the corresponding NLS profile. 




Figure 8: (color online) Comparison of the NLH solutions with and without the adjustment of incoming beam 
described in Section [6] (dashed red lines and solid blue lines, respectively), and the NLS solution (dotted black line). 



7.4- Comparison with the Previous Method 

In the nested iteration scheme of 0, 20, 21], at each outer iteration the Kerr nonlinearity is 
considered fixed, or frozen, which yields the linear homogeneous variable coefficient equation 

A + k^ + efcg|^(j)|2-) ^{J+i) = 0. (60) 

Equation (j60p is also solved iteratively, by building a sequence of Born approximations. In doing 
so, at each inner iteration an inhomogeneous linear constant coefficient equation 



(A + kl) = -ekllE^^'^'^l^'^E^^- 



■l,k) 



(61) 



is solved using the separation of variables. We will call this approach the "nested iterations method." 
The efficacy of this method can be improved by getting rid of the inner iterations ()6ip and solving 
equation (jGOp by the Gaussian elimination. We will call this the "freezing iterations method." 

In the one-dimensional case of [19|], the freezing iterations diverged above a certain nonlinearity 
threshold, while Newton's iterations converged for the entire range of nonlinearities of interest. In 
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the current multi-D cases that correspond to the critical NLS, i.e., D = 3, a = 1 and D = 2, a = 2, 
both the nested iterations method and the freezing iterations method diverge when the NLS solution 
collapses, i.e., when the input power is above Pq, while Newton's algorithm converges, at least for 
some configurations, thereby showing a much better efficacy. 

Another case of interest from the standpoint of applications is the subcritical NLS, D = 2, a = 1, 
which admits solutions in the form of spatial solitons. To compare the three methods in this case, we 
use each of them to repeat the simulation of Section [7.1.11 while varying the domain size ^'max- The 
quantity of interest is the threshold value Z^jnx = '^max'^^°^'^i below which a given solver converges 
and above which it diverges. The results are given in Tabled) We can see that the nested iterations 
method of 0, 0, ^ converges only for relatively short domains ^max < -^max'^'^"^'^ — ^2. Replacing 
the inner iteration by a direct solver brings along a certain improvement: ^max < ■^max'^'^"^'^ — 135. 
However, similarly to the one-dimensional case, Newton's iterations converge for the widest selection 
of cases, at least until Z^ax = 500. Moreover, this limit is due to the memory constraints rather 
than divergence, and the actual •^max'^^°^'^ may be even larger. 



Table 4: A comparison of the efficacy of tfie three methods for the soliton case D — 2, a = 1. Each method converges 
for Z^ax < ^mar"""' and diverges for Z^^x > Z'J^l^'"''""' . 



Method 


nested freezing (160p. (I6ip 


freezing (160p. LU solver 


Newton's 


^threshold 
rnax 


42 


135 


> 500 



8. Discussion and Future Plans 

In this study, we propose a novel numerical method for solving the scalar nonlinear Helmholtz 
equation, which governs the propagation of linearly polarized monochromatic light in Kerr di- 
electrics. The NLH is the simplest model in nonlinear optics that allows for the propagation of 
electromagnetic waves in all directions and, in particular, for backscattering, and accounts for 
nonparaxial effects. Our key result is that the NLH eliminates the singularity that characterizes 
solutions of the nonlinear Schrodinger equation, which is a reduced model based on the paraxial 
approximation. Another important finding is the discovery of narrow nonparaxial solitons and the 
development of numerical capability for simulating their collisions. 

Mathematically, the NLH is an elliptic equation, and must be solved as a nonlinear boundary- 
value problem. This presents additional difficulties for both analysis and computations compared 
to the traditional treatment based on the NLS. The latter has a predominant direction of prop- 
agation and requires a Cauchy problem. Physically, we consider the propagation of laser light in 
a layered medium with interfaces across which both the linear and nonlinear components of the 
refraction index may undergo jumps. The presence of material discontinuities necessitates setting 
the condition that the field and its first normal derivative be continuous at the interface. 

To solve the NLH numerically, we develop a fourth-order finite difference scheme for one, two, 
and three space dimensions (in the latter case we assume cylindrical symmetry). Finite differences 
are chosen over other possible approximation strategies because of their simplicity and ease of 
implementation. Indeed, the geometry of the problem enables a straightforward discretization on 
a uniform rectangular grid. On the other hand, having a high order scheme is important because 
it alleviates the point-per-wavelength constraint for large domains and also helps resolve the small- 
scale phenomenon of backscattering. In particular, high order accuracy must be maintained across 
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the material discontinuities. This is achieved by using special one-sided differences. In doing so, to 
simplify the overall discretization we move the outer boundaries away from the interfaces so that 
the artificial boundary conditions do not "interfere" with the interface treatment. The scheme used 
in the interior is of a semi-compact type, it is written on three nodes in the longitudinal direction 
and five nodes in the transverse direction. Having a compact three-node stencil in the longitudinal 
direction greatly simplifies both the treatment of the interfaces (no special "near interface" nodes) 
and the treatment of the outer boundaries (no non-physical evanescent modes). At the same 
time, a compact stencil in the transverse direction is not required because there are no material 
discontinuities. This circumstance greatly simplifies the design of the overall scheme. 

The second key component of the proposed algorithm is the nonlinear solver, which is based 
on Newton's method. The simulations of [l^ have demonstrated a clear superiority of Newton's 
method in the one-dimensional case. In this paper, we generalize our Newton's solver to the 
multi-dimensional case, with the expectation that it will let us solve the NLH for those settings 
when the NLS breaks down, namely, when the NLS solution becomes singular {cr{D — 1) = 2 with 
input powers above critical), or when the beam width becomes very narrow in the subcritical case 
(D = 2, cj = 1), or when counter-propagating nonpar axial solitons interact. 

The Newton's solver that we developed has indeed lived up to the promise. In the critical cases, 
it enables the central result of this work, which is the discovery of bounded NLH solutions for those 
cases when the corresponding NLS solution blows up. Physically, it shows that nonparaxiality can 
suppress the singularity formation and hence arrest the collapse of focusing nonlinear waves. While 
there may be other physical mechanisms that also help arrest the collapse (neglected along the way 
when the NLH was derived from the Maxwell's equations), it was not known until now whether the 
solution becomes regular already in the framework of the scalar NLH model, which is the simplest 
nonparaxial model that incorporates the backward traveling waves. 

Predictions of the NLH in the subcritical case include the existence of narrow nonparaxial soli- 
tons, and analysis of the interactions (collisions) of such beams, specifically in counter-propagation. 
These results may be of of relevance to potential applications, e.g., the design of the next generation 
of all-optical circuits. Note that in our previous work [i^] we have already been able to compute 
narrow spatial solitons. However, the new method proposed in this paper allows us to do that over 
much longer propagation distances, see Section mi 

Let us also note that a different configuration with counter-propagating solitons has been studied 
by Cohen et. al. in 22] using a system of coupled NLS equations which approximates the NLH. As, 
however, mentioned in 22], the coupled NLS model is not problem free as it is neither an initial- value 
problem nor a boundary-value problem. In contrast, since the NLH is solved as a boundary-value 
problem, it is a natural mathematical setting for such counter-propagating configurations. 

The computational cost of the proposed algorithm still remains relatively high; it is dominated 
(both in memory and CPU time) by the cost of inverting the Jacobian matrix using a direct method. 
This cost can be reduced if the LU decomposition is replaced with an iterative method. As, however, 
the Helmholtz operator subject to the radiation boundary conditions is not self-adjoint, the only 
viable choice of an iteration scheme will be a method of the Krylov subspace type. For this method 
to work, the system must be preconditioned, and it is the design of a good preconditioner that 
will be in the focus of our future work on the linear solver. Several candidate techniques will be 
investigated, including the constant coefficient Helmholtz operator to be inverted by the separation 
of variables and a paraxial preconditioner based on the Schrodinger operator. 

As far as the dependence of Newton's convergence on the domain size, see Section 17.2.1.21 we 
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attribute it to the generally known "fragility" and, in particular, sensitivity of Newton's convergence 
to the choice of the initial guess. On one hand, it is intuitively reasonable to expect that if the 
outer boundary is located in the region of maximum self-focusing, then the iterations may experience 
difficulties to converge, see Table [3l On the other hand, at the moment we do not have a clear 
and unambiguous mathematical explanation as to why exactly that happens. We have tried a few 
simple remedies, such as using a continuation approach in the nonlinearity coefficient e and using 
a damped NLS solution as the initial guess, but none of those has made a substantial difference. 



We note that in the one-dimensional case the exact solution was available in the closed form [19|] 
and hence we could at least test Newton's convergence by substituting this exact solution as the 
initial guess. In multi-D, however, we are not aware of any closed form solutions for the slab of 
finite thickness and therefore, a similar validation procedure becomes problematic. 

The piecewise constant formulation that we have considered in the paper in fact presents no 
loss of generality, at least from the standpoint of numerical solution. It can be very easily extended 
to the NLH with piecewise smooth material coefficients ^^{x) and e{x). All one needs to do is 
replace the constants ly and e in the definition of the scheme with the values at the corresponding 
grid nodes: i'n,m = i^izn, x±^rn) and en,m = ^{zn, x±^m)- However, while the resulting scheme will 
approximate the variable coefficient scalar NLH ([T]) with fourth-order accuracy, the validity of 
equation ([1]) itself from the standpoint of physics may be in question. Indeed, the derivation of 
the scalar NLH from Maxwell's equations in the case of variable coefficients introduces additional 
terms (spatial derivatives of v and e|-Ep) which are not included in equation ([T|). 

The layered structure and simple geometry that we have adopted present no substantial loss of 
generality, because this formulation corresponds to many actual physical (e.g., laboratory) settings. 
The plain-parallel setup studied in the paper certainly simplifies the discretization. At the same 
time, we are reasonably confident that the proposed scheme can be generalized to more elaborate 
geometries without compromising its high order accuracy, which is of key importance. One natural 
approach to doing that is to use Calderon's projections and the method of difference potentials 361 ]. 

From the standpoint of physics, the scalar NLH is certainly not the most comprehensive model. 
It is rather a reduced model based on a number of simplifications. Most notably, the vector nature of 
electromagnetic field is not taken into account by the scalar NLH because of the assumption of linear 
polarization. Vectorial effects, on the other hand, are known to become important close to when 
the nonparaxiality does, i.e., once the beam width becomes comparable to the carrier wavelength. 
Moreover, the scalar NLH governs monochromatic fields (continuous- wave laser), whereas the actual 
fields are always time-dependent (typically, pulses of certain duration). Nonetheless, if the duration 
of the pulse is sufficiently long (many oscillation periods), then the time-periodic model will provide 
a good approximation. 

To take into account the entire range of relevant physical phenomena one needs, of course, to 
go back and solve the full nonlinear Maxwell's equations. This, however, is a very challenging 
computational task and besides, the solutions of full Maxwell's equations may be hard to analyze 
or verify precisely because of all too many additional physical effects. That's why the analysis of 
the simplest nonparaxial model (i.e., the NLH) may provide a very useful insight into the relevant 
physics as, in particular, it allows to study the important phenomenon of nonlinear backscattering. 

Given the previous considerations, we believe that in the context of physics, the next most 
natural and most beneficial extension of the work presented in this paper will be taking into account 
the vectorial effects. The current work provides a solid foundation for this extension as many 
key elements of the algorithm, e.g., the nonlocal artificial boundary conditions, will only require 
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technical rather than conceptual changes. On the pure numerical side, in addition to the previously 
mentioned major modifications to the linear solver, we can consider a number of strategies aimed 
at further improving the numerical resolution in the regions of foremost interest (e.g., around the 
maximum self- focusing) while not increasing the overall computational cost. Examples include 
local grid requirement and/or combined approaches when most of the domain is to be done using 
the NLS whereas the local area of collision between the solitons is computed using the NLH. 



A. Continuity Conditions at Material Interfaces 



For optical frequencies, we can disregard all magnetization effects in the medium (see 37|, 
Chapter IX]) and write down the time-harmonic Maxwell's equations as follows: 

— B = cuvlE, -—D = cmlB, (62) 

c c 

where the specific form of how the electric induction D depends on the field E is not important for 
the derivation of the interface conditions. Note, however, that as our medium is a dielectric, both 
fields E and B, as well as the induction D, remain finite everywhere including the interfaces. 
Let an interface plane be normal to the coordinate z of the Cartesian system {x,y,z). Then, 

the first equation of (f62|) implies that the quantity {cutIE)^ = is bounded at the 

interface. As the derivative which is taken along the interface, is bounded in its own right, we 

conclude that -g^ is bounded. This immediately yields the continuity of Ey across the interface. 
The continuity of E^ can be established the same way, by taking into account the boundedness of 
(curl E)y = — Altogether, this means that the tangential component of the electric field 
E must remain continuous. Likewise, the continuity of the tangential component of B across the 
interface can be derived by employing the second equation of (f62]) and the boundedness of D. 
Next, consider the case of linear polarization: 

E = [E^, 0,0] and B = [0,By,0]. 

Then, the continuity of By immediately implies the continuity of because from the Faraday 
law (the first equation of (j62p ) we now have: ^By = Altogether, we conclude that for the 

linearly polarized light propagating through a (transparent) dielectric with material discontinuities, 
both the electric field E and its first normal derivative must be continuous at all the interfaces. 



B. Notation for Central Difference Operators 

We denote the central difference operators by the letter D with the order of differentiation in 
the subscript and the order of accuracy in the superscript. The full list for the finite differences in 
the X (or p) direction is as follows: 

^(2)^ dcf En,n.+l-^En,n.-l 
T^(2) dcf En 

^xxx^ = = OxxxEn,m + U [tl j. 
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n(2) Z? 42? En^rn+2 ~ 4-E'n,m+l + ^En^m ~ '^En^m-l + En^rn-2 _ ^ p _i_ ( h'^\ 
^xxxx^ — '^xxxxi^n,m + CI/ (^/l j, 

,m— 1 ,m— ^ 'in , /-/I /'il4\ 

-f^i -c/ = ■ — Y2h~, ~ (Jxi^n,m + ^ \ri j , 

(4-) def — £'n,m+2 + 16ii^n,m+l — SOE'n.m + IGE'n.m-l — ^n,m-2 r^ , ^/i/,4\ 

J->xxE = ^^^^2 = OxxEn,m + \ri ) . 

Because of the semi-compact approximation we use, only the second-order operator is required 
in the z direction 
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